TIME INTEGRATION AND STEADY-STATE CONTINUATION FOR 
2D LUBRICATION EQUATIONS 

PHILIPPE BELTRAMEtt AND UWE THIELEt§ 

Abstract. Lubrication equations allow to describe many structurin processes of thin liquid 
films. We develop and apply numerical tools suitable for their analysis employing a dynamical 
systems approach. In particular, we present a time integration algorithm based on exponential 
propagation and an algorithm for steady-state continuation. In both algorithms a Cayley transform 
is employed to overcome numerical problems resulting from scale separation in space and time. An 
adaptive time-step allows to study the dynamics close to hetero- or homoclinic connections. The 
developed framework is employed on the one hand to analyse different phases of the dewetting of a 
liquid film on a horizontal homogeneous substrate. On the other hand, we consider the depinning 
- - ■ of drops pinned by a wettability defect. Time-stepping and path-following are used in both cases 

' to analyse steady-state solutions and their bifurcations as well as dynamic processes on short and 

' long time-scales. Both examples are treated for two- and three-dimensional physical settings and 

, prove that the developed algorithms are reliable and efficient for Id and 2d lubrication equations, 

04 1 respectively. 
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I— i' 1. Introduction. The dynamics of structuring processes of thin Hquid films, 

ridges and drops on solid substrates can be described by so-called thin film or lubri- 
I cation equations. The notion 'thin' refers to systems where the mean thickness of the 

■ film/drop is small as compared to all typical length scales parallel to the substrate. 
[ This allows for the usage of the long wave approximation |65] . Different varieties of 

, ^ such an equation may model, for instance, dewetting due to van der Waals forces 
I [751 [Ml [8il 1101^ [i]. the evolution due to a long- wave Marangoni instability of a film 

■ heated from below [66l [l2l ES] , and the instability of a film of dielectric liquid in a 
I capacitor [56^ I108[ |46] . Including driving forces parallel to the substrate allows 

to describe, e.g., droplets that slide down an incline under gravity under isothermal 
I [69l I98j and non-isothermal conditions [El [95], the evolution of transversal front in- 

■ stabilities [571 [3 [Ml [13 [M], and the evolution of shocks in films driven by a surface 
I tension gradient against gravity [101 [88]. A large number of extensions exist to de- 

■ scribe, for instance, two-layer films, films with soluble or non-soluble surfactants, films 
I of colloidal solutions or to incorporate effects of evaporation, complex rheology or slip 

■ at the substrate. For reviews see, e.g. references |65 [ [89 l [49 l [99 l I13j. 
^ I Thin film equations are extremely versatile in their application and are related to 

other 'standard' equations employed in studies of pattern formation out of equilibrium 
. , [22]. It can, for instance, be shown that the (convective) Cahn-Hilliard and the 

■ Kuramoto-Sivashinsky equations may be obtained from a thin film equation for sliding 
drops and fiowing films as limiting cases for small and large lateral driving, respectively 

m- 

Due to the strong non-linearity a typical thin film equation is difficult to handle 
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numerically, in particular, if the physical situation that shall be described is three- 
dimensional resulting in a partial differential equation (PDE) with two spatial dimen- 
sions (2d). In the following we will exclusively employ the latter counting style of 
dimensions, i.e., a two-dimensional physical situation resulting in a PDE with one 
spatial dimension will be referred to as "one-dimensional case (Id)". 

The present work develops and applies a common numerical framework that can 
be used to study the time evolution and to follow steady state solutions in parameter 
space in both - the Id and 2d - cases. For viscous fluids (small Reynolds number 
flows), for a surface tension that dominates over viscosity (small capillary number), 
and a small lateral driving force (e.g., small inclination angle of an incline or small 
temperature gradient) the long wave approximation results in the following equation 
for the time-evolution of the film thickness or droplet profile h{x,y,t) |65[ l49] 

(1.1) dth = -V • [m{h)\Ip{h) + ^l{h)] , 

where m[h) is a mobility function, fl{h) represents the lateral driving force, and the 
pressure p{h) may contain several terms that determine the dynamics. A curvature 
or Laplace pressure results from capillary and stabilizes a flat fllm. Its contribution 
results in a Bilaplacian of the height h in Eq. I|l.ip . Because it is the highest or- 
der derivative in the equation it constitutes one of the main numerical difficulties. 
Pressure contributions that destabilize the flat fllm can result from various physical 
mechanisms [65l [25] . Examples include electrostatic pressure terms for dielectric Hq- 
uids in a capacitor [56l [60l fT08] . a pressure resulting from a magnetic fleld acting on 
a ferromagnetic fluid [H [14], a disjoining or conjoining pressure for very thin fllms 
below 100 nanometers thickness resulting from effective molecular interactions be- 
tween the substrate and the free surface (wettability effects) [24} W!\ [43] , a 'thermal 
pressure' for a thin fllm on a heated plate (where the long-wave Marangoni mode is 
dominant) [106[ [6il [T2] and a hydrostatic pressure due to gravity, e.g., for a fluid 
fllm under a ceiHng [32l [26l [16] . All the mentioned destabilizing effects result in a 
long-wave instability, i.e., an instability with wavenumber zero at onset; similar to the 
Rayleigh- Taylor instability. 

In the present work we use two variants of a disjoining pressure as generic exam- 
ples for a destabilizing mechanism because their particular thickness dependencies are 
numerically demanding and extensive literature results allow for detailed comparison. 
Note, however, that the developed algorithms for time integration and continuation 
can be readily appHed to any other combination of stabilizing and destabifizing pres- 
sure terms. Using a disjoining pressure formulation allows to consider the presence of 
an ultra-thin wetting layer on the macroscopically 'dry' parts of the substrate. Such 
a fllm is often called a 'precursor fllm' [2l]. The existence of such a fllm removes the 
contact fine singularity [3Tlll02| . 

In the absence of a lateral driving force (^ = 0) , a fllm that is linearly unstable 
evolves during a short-time evolution into a structure of holes, drops or mazes with 
a typical structure length determined by fllm thickness and other control parameters 
that reflect the character of the destabilizing phenomenon [73 l [84 l [64 l [52 l [T H [H j [81] . 
This process is often called 'spinodal dewetting' [61]. However, the resulting short- 
time structure is unstable (representing a saddle in function space) with respect to 
coarsening and on a large time-scale the structures coarsen until the system eventu- 
ally approaches the global energetic minimum, i.e., a single drop or hole [8l [36t [91]. 
Indeed without a lateral driving force the system follows a relaxation dynamics and 
the evolution equation can be written using the variation of an underlying energy 
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functional [Ml EH [HU [99] . The variational form of the (then coupled) evolution equa- 
tions can also be derived for multilayer films in similar settings [70l[7l]. However, 
even in the one-layer case, details of the different phases of the process are still un- 
der investigation. Examples include the initial structuring process that might occur 
via nucleation or a spinodal surface instability [101\ l4]. and the mechanisms of the 
transverse instability dewetting fronts [85tl7i]. 

A detailed understanding of dewetting films and the other processes listed above 
is only possible if the pathways of time-evolution and the steady states described by 
Eq. I|l.ip can be determined in the one- and two-dimensional case using a fast and 
versatile algorithm. As the steady droplet solutions in Id can be seen as periodic 
trajectories in a conservative dynamical system one can use available continuation 
packages for ordinary differential equations [28^ \29\ [30] to map different solution fam- 
ilies and their linear stability (see, e.g., [103[ [95 l llOOj ). Careful interpretation allows, 
e.g., to predict for dewetting films the dominance of different rupture mechanisms 
within the linear unstable parameter range |101l I103[ [89] . No such continuation tools 
are, however, available in the two-dimensional case. 

The situation is more involved when lateral driving forces are present, i.e., gravity 
on an incline [i2l 11071 [69] . or temperature gradients along the substrate [TSl [50l [9]. 
There, new phenomena appear like transversal instabilities of sliding liquid ridges, 
advancing and receding fronts |42[ [87j \7\ [SH [98l [94] . Another fascinating finding is 
related to sliding drops: Beyond a critical driving force the drop forms a cusp at the 
back end and 'emits' smaller satellite droplets [69t [5} [54]. For a driven contact line, 
heterogeneities of the substrate can cause a stick-slip motion [80] . In the setting of a 
sliding droplet this leads at a critical driving force to the depinning of droplets from 
such localized heterogeneities. In the vicinity of the responsible sniper bifurcation 
the resulting motion resembles stick-slip motion: The drop sticks a long time at a 
wettability defect and then suddenly slips to the next defect. This was studied in the 
Id case in Refs. [93 [96]. Differences in time scales for the stick- and the slip-phase may 
be many orders of magnitude. On the homogeneous substrate one can, in the Id case, 
regard stationary periodic droplet and surface wave solutions as periodic trajectories 
of a dissipative dynamical system that emerge from the trivial fiat film state via a Hopf 
bifurcation [93- This allows to employ standard continuation packages |30j to obtain 
solution families and to track the various occurring bifurcations [102[ [45] I100[ [92] . 
The same applies for fronts (shocks) and droplets that correspond to heteroclinic and 
homoclinic orbits, respectively [TOl [Ml [H] • Note, however, that at present little is 
known about the solution and bifurcation structure in the 2d case. 

Having all the described physical applications in mind, we conclude that for set- 
tings with and without lateral driving forces a powerful algorithm that allows to follow 
steady state solutions in parameter space and to simulate the time evolution would 
be highly beneficial. Here, we seek to develop and to apply such an algorithm. In 
particular, we aim at a time integration scheme with an adaptive time-step and tools 
for bifurcation analysis that are applicable in the Id and (most importantly) 2d case. 

Semi- implicit time-stepping schemes (cf. Ref. |109j ) are a standard method to 
perform the time- integration of non-linear equations. The method treats a linear 
operator implicitly and the remaining non-linear terms explicitly. The method is 
successful, e.g., when dealing with the Navier-Stokes equations |105) . because the 
fastest time scale normally arises from the linear Laplacian (or Stokes) operator. In 
the case of the lubrication equation l|l.ip , the natural candidate for the linear operator 
seems to be the bilaplacian. Indeed, it is very relevant for small perturbations of a 
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flat film. However, for states like a drop that coexists with an ultra-thin precursor 
film scale-separation makes this approach problematic. Alternatively, to avoid the 
problem one may choose to take at each time step the Jacobian matrix as the linear 
part. In consequence, the semi-implicit scheme requires to solve a linear system whose 
conditioning depends on the one of the Jacobian matrix. For an equidistant spatial 
discretization, we expect that the Jacobian is very badly conditioned (worse than in 
the case of the Stokes operator) because of the presence of the fourth order bilaplacian. 
Although, specialized algorithms exist that economically solve the linear system when 
the Laplacian is the implicit linear operator, no such algorithm seems to exist for 
general asymmetric Jacobian matrices. Another way that avoids this problem might 
be the use of a (pseudo-)spectral method [HlIIl]- However, the presence of higher 
derivatives contributes to a large spread in the eigenvalue spectrum which may lead to 
numerical problems. This effect becomes even larger because of the particular solution 
types we want to consider: Sliding drops on homogeneous or heterogeneous substrates, 
in particular, the cusp formation [69] lead to a further enlarged eigenvalue spectrum. 
Therefore we do here not pursue such a way. Instead we focus on an algorithm 
based on a spatial finite difference discretization on a uniform mesh of N nodes. The 
resulting large and sparse discrete system of equations contains obviously strong non- 
linearities. We seek an eflective iterative method to carry out the integration in time 
and the tracking of steady states in parameter space. 

An efficient time integrator can be based on an exponential propagation scheme 
that does not require any preconditioning and is known to be stable even for stiS non- 
linear systems. The scheme is based on the exact solution of the linearized equation at 
each time-step requiring the computation of the exponential of the Jacobian matrix. 
This operator is not directly computed but, as commonly practiced for large and 
sparse matrices, only its action on vectors is estimated using projections on small 
Krylov subspaces of dimension K <^ N [76j. The standard algorithm to perform this 
task is the Arnoldi procedure possibly incorporating improvements as proposed by 
Saad [77] ■ As a result the exponentiation can be performed at a negligible cost as 
long as K is not too large. 

The application of such a scheme to lubrication equations proves to be reliable. 
However, the approximation in Krylov subspaces converges rather slowly. We show 
that the necessary dimension K is about one hundred while in Refs. [34l [77l l40j I104j 
fiT « 10 is sufficient. The difference results from the presence of the fourth order 
Bilaplacian in lubrication equations. Its effect is more disadvantageous than the one 
of a second order operator because the magnitude of the large negative eigenvalues 
increases with the order of the differentiation operator [55]. To our knowledge, no 
numerical study of the efficiency of the Krylov subspace approximation to a matrix 
exponential operator that contains a fourth order operator exists in the literature. 
Here we propose to improve this step by coupling the already developed methods 
to determine the rightmost spectrum with the classical Arnoldi procedure [59]. In 
particular, the application of a Cayley transform proves to be the most powerful 
method. The resulting scheme is particularly suited to efficiently adapt the time-step 
to the corresponding time-scale of the dynamics. This allows for very large time-steps 
as, for instance, in problems involving coarsening and stick-slip drop motion. 

The outlined scheme can not only be used for time-stepping. We show that the 
Krylov reduction associated with the Cayley transform can as well be applied to track 
steady states in parameter space employing a continuation scheme that consists of the 
determination of a tangent predictor, the application of Newton's algorithm along the 

4 



secant direction, the detection of bifurcation points and finally the determination of 
the direction of bifurcating branches of steady solutions [82] . 

The paper is structured as follows. Section [2] presents the lubrication equation 
and its spatial discretization. Then we describe in Section [3] the exponential prop- 
agation method and discuss different ways to exponentiate the Jacobian before in 
Section [H we adapt the developed algorithms to employ them for the continuation of 
steady states. Two appendices give details for the Krylov reductions Q and the im- 
plementation and convergence of the algorithms (|B]). In the remaining part, we apply 
the developed numerical algorithms to two typical physical situations. The first one 
(Section [5]) is the dewetting process of a thin film on a horizontal solid substrate. The 
slow coarsening process that follows the initial fast patterning provides an excellent 
test for the adaption of the time-step. We do as well study steady state solutions in 
2d, in particular, we discuss several solution branches corresponding to quadratic and 
hexagonal arrays of drops. Second, in Section[6]we study the pinning/depinning prob- 
lem for ridges and drops on an inclined heterogeneous substrate, where a wettability 
defect can retain a liquid ridge or drop up to a critical inclination. The path-following 
algorithm is appHed to determine branches of steady state solutions and, in conse- 
quence, the onset of depinning in the 2d case. The stick-slip motion of drops beyond 
depinning is investigated utilizing our time-stepping algorithm. For comparison with 
the literature we do as well provide selected results for the Id case. Our conclusions 
are found in Section \7\ 

2. Modelling and Spatial Discretization. 

2.1. Lubrication equation. Consider a liquid layer on an (inhomogeneous) 
one- or two-dimensional solid substrate, corresponding to a two- or three-dimensional 
physical situation, respectively (Fig. ??). The Hquid partially wets the substrate and 
might be subject to a constant lateral force P. Using the long-wave approximation, 
the dimensionless evolution equation for the film thickness profile h{x,y,t) derived 
from the Navier-Stokes equations, continuity and boundary conditions reads [65^ l49] 

(2.1) dth ^ -V • {m{h) [V {Ah + n{h, x)) + Pe^]} , 

where V = {dx, dy) is the planar gradient operator and A = d^^ + dyy is the planar 
Laplacian. Note, that Eq. I|2.ip has one and two spatial dimensions for two- and 
three-dimensional physical settings, respectively. In the following they are referred to 
as the Id and 2d case, respectively. The mobility function m(/i) = corresponds to 
Poiseuille fiow without slip at the substrate. The term A/i represents the Hnearized 
Laplace pressure for small slopes of the film thickness profile. The wetting properties 
are modeled by the disjoining pressure Ii{h^x). Thus, for a striped heterogeneous 
substrate the disjoining pressure will depend not only on film thickness but also on 
the .T-direction. The lateral force acts as well in the x-direction (Fig. ??). Many 
particular forms of the disjoining pressure are discussed in the literature. The most 
common ones that are used in simulations model the presence of an ultra-thin wetting 
layer (a so-called 'precursor film') of about 1-10 nm thickness. Note that, the resulting 
short-time dewetting dynamics is often called initial 'film rupture'. However, the 
disjoining pressure allows for a stable precursor film there is no 'true' film rupture. 
The short-time dynamics results in large amplitude structures, either drops coexisting 
with a thin precursor film or holes in the thick film that do, however, not reach down 
to the dry substrate but are still covered by the precursor film. We will employ the 
notion 'rupture' in this meaning at several points in the paper. 
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To facilitate comparison to the literature here we employ the two expressions for 
the disjoining pressure: 



(2.2) U{h) = -2e-''{l-e-'')-Gh 
(as used in [TQ2tfT03]). and 

(2.3) u{h) = -^ + e-'^ 

(as used in [83l I101| ). We call them (I) and (II), respectively. Varying wettability 
properties due to a heterogeneous coating can be introduced in the short-range part 
of the disjoining pressure. We employ case (II) and have 

(2.4) uih,x)^^-[l + eax)]e-\ 

where ^{x) is the heterogeneity profile and e the amplitude of the heterogeneity (as 
in [93 ES]). 

2.2. Functional space. The considered domain is D = [0', L] (Id case) or D = 
[0;Lx] X [0;ij,] (2d case). Eq. (|2.ip defines a partial differential equation (PDE) in 

(2.5) dth = F{h,x),h e H^iD) 

with F{h) = -V-{m(/i) [V {Ah + U{h)) + Pe^]}. In order to obtain a well posed PDE 
system, we introduce periodic boundary conditions. In the studied case of non- volatile 
liquids the mass M = h{x,y)dxdy is conserved. The measure H — M/ J^dxdy 
represents the mean height and u = h — H is the perturbation. The PDE (|2.5p defined 
in the space 

Eo^l^u e H^{D) : J u{x,y)dxdy = 

(2.6) and periodicity of and 0<i,j<3} 

is well defined: (i) F operates on the EucHdian space {H + Eq}, and (ii) linear 
operators operate on the Hnear space Eq. 

Note that we do not impose a positivity condition {h > 0) in the functional 
space. Indeed, Zhornitskaya and Bertozzi present a positivity-preserving nu- 

merical scheme for lubrication equations which do not contain a disjoining pressure 
(n = 0). The result is based on the fundamental result of the existence of entropy 
preserving solutions [6]. In the case of the presence of a disjoining pressure, the posi- 
tivity of such schemes is guaranteed only for the Id case [HIS]. Nevertheless, entropy 
consistent finite-element schemes can be developed for the 2d case [38]. Here, the 
presence of a precursor film makes the positivity preserving problem less crucial. Fur- 
thermore, the error estimator developed with our scheme is sufficiently accurate to 
prevent negative solutions or 'numerical' singularities. Note finally, that the concept 
of entropy solutions might not be applicable to films exposed to a driving force. 

2.3. Spatial discretization. A finite difference scheme is used associated with 
a regular meshing of the substrate domain. We use in x and y directions, (N^ + 1) 
and {Ny + 1) mesh points at distances Sx = N^/Lx and Sy — Ny/Ly, respectively. 
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The number of discretization points is denoted hy N ~ x Ny in the 2d case and 
N = Nx + I in the Id case. The differentiation operators are approximated by a 
standard centered five-point stencil. This means that for the Bilaplacian and the 
third order operator VA it is a second-order approximation while for the gradient 
or divergence operator V and the Laplacian A, it is a fourth-order approximation. 
Thus, the differentiations operators are sparse band matrices with maximally 5 (Id) 
or 25 (2d), non-zero elements on each fine. The map F{h,x) defined in Eq. I|2.5p is 
discretized according to Eq. I|2.ip . i.e., in the mass conservation form. Our algorithm 
involves as well the Jacobian of the map F, which is discretized using Eq. I|2.8p . This 
results in a second-order approximation. Thus, the discretized Jacobian matrix is not 
in a conservative form. With this choice we avoid a further increase of the number 
of non-zero elements. Indeed, in the conservative form the matrix associated with 
the divergence operator would multiply the other differentiation matrices. However, 
the use of a matrix in a non-conservative form increases the mass fiuctuation of the 
vector h due to numerical round off errors. This may result in a significant variation 
of the mass after several time-steps. A solution could be to integrate the restriction 
of the E space directly into the operator. Then, however, one would loose the band- 
structure of the operator that constitutes a property of paramount importance for 
sparse matrix manipulations like, e.g., the incomplete LU-factorization. Therefore, 
we keep the matrix without additional modifications and impose mass conservation 
during the Krylov reduction as explained in Section [Aj 

2.4. Jacobian matrix. Next, we justify our choice to compute the Jacobian 
matrix at each time-step. First, we explicitly determine the Jacobian matrix as lin- 
earization of the function F{h, P) at /i = ho 

(2.7) 3 = DhF{ho,P) 

where Dh denotes the differentiation operator with respect to h. The Jacobian matrix 
is a sum of differentiation operators up to fourth order: 

(2.8) Ju = VQ{hoj x)u + vi{hQ,x) ■ 

-t-moIIoAu — Vmo • VAm — moA^u. 

with 

(2.9) vo{hQ,x) = -moA^/io - Pd^rn'o 

+moAno + Vmg • VEq + Vm • VEg + mAU'o 

(2.10) vx{ho,x) = -m'^VAho + m'f^VUo + 2moVn' + VmoW + Pe^ 

and 

mo = m{ho) m'o = ^{ho) 
Ho = n(/io,x) U'o = ^{ho,x) 

The Jacobian matrix J has eigenvalues with large negative real parts. They correspond 
to the fastest time scales related to the presence of the differentiation operators. As 
these eigenvalues are situated in the left half-plane of the complex plane they are 
called "leftmost eigenvalues". The presence of such eigenvalues is the main cause 
for numerical spatial oscillations occurring when explicit methods are used for time 
integration. The resulting stability restriction on the time-step is overcome by treating 
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the linear term implicitly. Then, the main difficulty is to solve the resulting linear 
system. The leftmost eigenvalues are due to the differentiation operator of the highest 
degree. Therefore, on can filter out the leftmost eigenvalues by treating only this 
operator implicitly. In the literature the implicit linear term normally corresponds to 
the Laplacian for second order equations like, e.g., reaction-diffusion or Navier-Stokes 
equations and to the bilaplacian for fourth order equations Hke, e.g., the Kuramoto- 
Sivashinsky equation [35]. In our case, a natural candidate is the bilaplacian. However, 
it is multiplied by the Poiseuille fiow mobility mo = Hq. This implies that, e.g., for 
a height variation by a factor 10 like for a drop solution, a factor 10'^ appears in the 
linear operator. The resulting spatial scale separation complicates the situation as 
compared, e.g., to the Kuramoto-Sivashinsky equation. 

Furthermore, here the other differentiation operators may contribute to the left- 
most eigenvalues because the vectorial factor in front of these operators can have 
very large elements. For example, the Laplacian operator is scaled by tooIIo what is 
related to the disjoining pressure that can be very large close to the edge or contact 
region of drops. In this situations the Laplacian operator results in a non-neghgible 
contribution to the leftmost spectrum of J. 

In conclusion, we find that it is not advisable to filter the leftmost eigenvalues 
by employing a constant linear operator. It is preferable to compute the Jacobian 
matrix through a linearization at each time-step. However, the resulting asymmetric 
non-constant Jacobian matrix is not the best choice in an implicit integrator: The 
linear system is difficult to solve as the Jacobian is badly conditioned and there exists 
no general preconditioner method to improve the iterative methods. An alternative 
way is to employ exponential propagation as then there is no need to solve a linear 
system. Furthermore, it is known in the literature that exponential propagation is 
well suited for time evolution problems involving strong non-linearities. 

3. Time integration. 

3.1. Exponential propagation and error estimator. Starting from the known 
profile ho at to , the exponential propagation scheme consists in solving the autonomous 
evolution problem 

(3.1) ^^Fih,x), 

(3.2) h{to) = ho 

at each time-step. The approach for solving this system is to develop the operator 
F{h, x) near the state ho in a Taylor series 

(3.3) F{ho + u,x)= F{ho, x) + DhF{ho, x)u + R{u), 

where DhF(hQ,x) is the Jacobian matrix at the point ho and R{u) = O(it^) contains 
the quadratic terms. In order to simplify the notation we let 

(3.4) 3 = DhF[ho,x) 

(3.5) b^F{ho,x). 

The height variation u = h — ho is then the solution of the evolution problem 

fill 

(3.6) — =b + Ju + R{u), 

uito) = 
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which admits 

(3.7) u{T)^G{JT)bT+ f cxp((T-r')J)i?(M(r'))dT', 

Jo 

with 

(3.8) G(Jr) = £^^pM^ 

Jr 

as solution. The last expression is formally correct even if J is not invertible since 
the operator G can be expressed as the series G{x) = X]n=^-^"/(" + ^V-- Since 
the functions G and exp tend to zero at — oo, they filter the leftmost spectrum of 
the operator J. For a semi-implicit scheme one has a similar situation: There is a 
rational function of J that plays the role of the filter. Yet, a rational function tends 
slower to zero than G and exp when approaching — oo. Therefore, we expect the 
exponential propagation scheme to better filter the leftmost eigenvalues than the semi- 
implicit scheme. We will again take up the discussion of the advantage of choosing 
an exponentiation scheme over a semi-implicit scheme in Section [BJ Note, that the 
exponential propagation scheme has the important advantage of the existence of an 
analytical error estimator. Differentiating p.7p one obtains 

(3.9) dtu{T) = cxp(Jr)6 + [ cxp((T - r')J)dr' + i?(M(T)) 

Jo 

allowing to deduce the analytical error 

u, = T[dtiu){t) - F{h + u)], 
and to define the relative error 

(3.10) er = -. 



h + h\\' 

The latter can be used to efficiently control the numerical error. Close to equilibrium, 
however, the variation u is small and one has to introduce an absolute tolerance 
The error criterion becomes 

' l + ea\\u + h\\ 

which is similar to criteria proposed by [31] and |109j . 

3.2. Order of the scheme. If we consider only linear variations of the operator 
F{h,x), i.e., R{u) = in Eq. p.3p . the first term of Eq. I|3.7p is the exact solution of 
the evolution problem 13.61 

(3.12) u(i)(r) = G(Jr)6T. 

The solution h^^^r) = u'^-^^t) + ho is then an approximation of second order with 
respect to time. If one wants to employ an exponential schemes of a higher order one 
has to take into account the non- linear term R{u) in Eq. I|3.7p . Then the perturbation 
u can not any longer be obtained explicitly from Eq. (|3.7p . The strategy developed 
in [34] is to estimate u by successive approximations of the non-linear terms through 
the series u^^\ Then, the resulting scheme of order ^ -I- 1 is given by 

(3.13) u(^) (r) = (r) + J] cxp((T - r')J)4 (-) dr' , 
where the computation of the vectors c^. is detailed in [Sj . 



3.3. Krylov projection. Since Ref. [3l] was published a plethora of high order 
approximations were developed [77 l l40l |4T1 1104| . However, the accuracy of the expo- 
nential scheme does not only depend on the non-linear approximation of R{u) but as 
well on the approximation of the vectors Vg = G(Jr)& and = cxp(Jr)c. 

The latter are commonly performed using a projection onto a small Krylov sub- 
space of dimension m computed by the classical Arnoldi algorithm. In the literature, 
this step does not constitute a difficulty as a good approximation of the action of the 
above operators is already obtained with a small Krylov subspace of about m ~ 10. 
In contrast, in our case we find as detailed in appendix[Al that when using the Arnoldi 
method to obtain accurate results the Krylov subspace has to have a dimension of 
about TO = 100. How can the bad convergence in the case of a thin film equation be 
explained? All examples chosen in the literature are second order equations (reaction- 
diffusion equation, Schrodinger equation). The thin film or lubrication equation, how- 
ever, contains a fourth order operator. This indicates that the bad convergence of the 
Krylov Arnoldi algorithm results mainly from the presence of the bilaplacian as it has 
a worse conditioning than the Laplacian. 

Therefore, one has to improve the Krylov projection step in order to be able to 
apply the exponential propagation method efficiently to lubrication equations. As our 
improvement can be applied in the same way to schemes of any order, we focus on 
the simplest exponential scheme, i.e. the linear approximation of order 2. 

We refer the reader to appendix [A] for a detailed description of the two methods 
used to estimate G(Jt)6 and exp(Jr)c - the classical Krylov- Arnoldi projection and 
our proposed variations of that method. The key-idea of our improvement is to 
transform the spectrum of the operator J in order to accelerate the convergence of the 
Krylov approximation. It is well known that the Krylov- Arnoldi algorithm first tends 
to the part of the spectrum that has the largest modulus. However, the rightmost 
eigenvalues of J are the ones of primary interest for the time-stepping. To reach fast 
convergence we need to apply a transformation that allows these "wanted" eigenvalues 
to become the ones of largest modulus. Such transformations are commonly used 
when estimating the rightmost spectrum (see e.g. [llOj ). One can distinguish two 
main methods: Chebyshev acceleration and Cayley transform [59]. Only the latter 
one proves to be efficient in our case (SectionEI- However, the Cayley transformation 
requires an Incomplete-LU (ILU) factorization which needs 0(iV'^/^) steps. It turns 
out that a Cayley-Krylov method is the most efficient one as long as the system size 
stays below w O(IO^). This is shown in the appendix. 

To summarize, the time integration of Eq. (|2.ip is performed using an expo- 
nential propagation scheme that employs Krylov projection. The scheme is stable 
independent of the particular Krylov approximation used. For moderate system sizes 
of A^ = O(IO^) the Cayley-Krylov projection furthermore allows to employ adaptive 
time-stepping. As it provides as well the leading modes at each time-step it proves to 
be a good tool to better understand the film dynamics. 

4. Continuation of steady-state solutions. 

4.1. Introduction. In this section we develop an algorithm to follow branches 
of steady-state solutions when varying a parameter p, i.e., we seek the branch {h,p) 
such that 

(4.1) F{h,p)^0. 

We adopt the 'classical' tangent predictor - secant corrector scheme [82] • For both 
steps the Jacobian has to be inverted. This operation is quite different from the ex- 
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ponentiation required in the time-stepping algorithm developed above. However, the 
procedure of Krylov reduction of the Jacobian using the Cayley transform is relevant 
for the needed inversion. Indeed, only a few eigen-directions associated to the eigenval- 
ues around the origin of the complex plane are relevant for the inversion. Eigenvalues 
of large modulus -namely, our leftmost eigenvalues- are of negligible influence. Fur- 
thermore, their presence can even lead to numerical oscillations. Therefore, the Cayley 
transform is a good algorithm for the determination of the wanted eigen-directions 
using a Krylov reduction. The resulting detailed knowledge of the leading eigenvalue 
facilitates the stability analysis and allows to detect bifurcation points. 

4.2. Tangent predictor - secant corrector method. First, we describe a 
continuation step as sketched in Fig. ??. One starts at point (/io,-Po) that represents 
a steady state ho associated to the parameter value Pg ■ Differentiating equation (|4.ip 
one obtains the tangent direction [ut,pt) of the continuing branch at the point (/iq, Pa) 
as solution of the system 

(4.2) 3oUt^~DpF{ho,Po)pt 

where Jo ~ DhFQiQ, Pq) is the Jacobian. The {ut,pt) solution is entirely determined 
by flxing the amplitude of pt . This is done by flnding the maximal amplitude of pt 
such that 

(4.3) \\Fiha+ut,Po+pt)\\<et\\ho\\. 

In this way one is not too far away from the steady state branch. Typically we take 
10""^ < Ct < 10~^ depending of the kind of problem. We denote this intermediate 
point by {hi, Pi) with 

(4.4) hi = ho+ut 

(4.5) Pi=Po+Pf 

The sign of pt remains to be chosen. It only changes at saddle-node bifurcations. 
In the (P, /(P)) plane the passing of a saddle-node bifurcation is characterized by a 
sign change of /'(P) and an increase of /'(P) before reaching the bifurcation. If both 
conditions are fulfilled the sign of pt has to be changed. 

Next, one uses Newton's method to solve Eq. I|4.ip close to (hi, Pi). The secant 
direction is taken orthogonal to the tangent direction {ut,pt) to ensure that one does 
not loose the branch of steady-state solutions even in the neighborhood of a turning 
point. For one of the Newton iteration steps that starts from {hk,Pk) the condition 



reads 

(4.6) JkWfc+i + DpF{hk, Pk)Pk+i = -F{hk, Pk) 

(4.7) {ut,Uk+i)=0 

(4.8) hk+i = hk + Uk+i 

(4.9) Pk+i=Pk+Pk+i 

where Jk — DiiF{hk, Pk) and Uk+i,Pk+i are the unknowns at fc + 1 Newton step. 



The Newton equation l|4.6p and the orthogonality condition l|4.7p can be written as 
the matrix equation 



Jk DpFiuk,Pk) ' 




Uk+l 




■ -F{uk,Pk) 


Ut Pt 












-F{uk,Pk) 
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(4.10) 



N 



Uk+l 

Pk+i 



One clearly sees, that the continuation step requires inversions of the Jacobian matrix 
J (tangent predictor) and of the matrix N (Newton corrector steps). Except at 
bifurcation points, the tangent equation (|4.2p is invertible in the space Eg. However, 
as explained before, we do not include mass conservation and continuity requirements 
of the space Eq into the Jacobian matrix in order to preserve its band structure 
because of its importance for a fast ILU-factorization. The restriction to the space 
Eq of the Hnear systems is ensured during the Krylov reduction (see Section EJ. As 
described next, here we apply the Cayley transform to perform the inversions. 

4.3. Estimation of the tangent predictor. The Cayley-Krylov reduction of 
the equation (|4.2p leads to 

(4.11) KiJmV;*,ut = -bpt = -DpF{ho, Po)pt. 

where Vm is the Krylov basis consisting of m vectors in Eq. It is constructed by 
letting the operator C = (J — cl)~^ act on the vector b = DpF{uo,po). The choice 
of the scalar c follows the same rules as in the time-stepping (section E}. Using the 
QR-method one now obtains the spectrum of the reduced m x m .Jacobian: 

(4.12) J^-ri = PjnZ^^^Pj^ . 

For the inversion we distinguish two cases: (i) if the kernel contains a non-zero eigen- 
vector vq, then the pair (uq, 0) is the solution of the problem; (ii) otherwise we perform 
the inversion and the solution is given by 

(4.13) ut = -KnPmD^-iP™"V„\6pt. 

In this way one is on the one hand able to detect bifurcation points and on the other 
hand one ensures that the continuation algorithm does not break down at turning 
points. 

4.4. Computation of secant corrector. The same method is applied to per- 
form one Newton step (|4.10p . However, the matrix DpF{hu, Pk) is not a band matrix 
because of the presence of the full vectors: Dp in l|4.6p and ut in the orthogonality 
condition (|4.7p . Since only J is a band matrix we solve the equation 

(4.14) Juk+i = -F{uk,Pk) - DpF{uk,Pk)pk+i. 
By solving the equations 

(4.15) 3usi = -F{uk,Pk), 

(4.16) 3us2 = -DpF{uk,Pk), 

the solution space of l|4.14p is the space of solutions of the linear combination of l|4.15p 
and l(4lll) 

(4.17) Mfe+l = Usl + Pk+lUs2- 

Finally, the scalar Pk+i is obtained by applying the orthogonality condition (|4.7p to 

Uk- 

(4.18) (lit, Usl) + Pk+l {ut, Us2) = 

Thus a Newton step needs an ILU-factorization of the Jacobian matrix and the solu- 
tion of two equations using the Cayley-Krylov method. The latter has negligible cost 
compared to the ILU-factorization. 
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4.5. Bifurcation analysis and stability analysis. Because the rightmost 
spectrum of the Jacobian is known we have all information needed to assess the 
stability of solutions. Furthermore, during the tangent predictor step one is able to 
detect the possible presence of a bifurcation. The direction of the bifurcating branch 
may be deducted from the eigen-directions of the kernel. Although, the rightmost 
spectrum is normally well evaluated it may occur that the rightmost eigenvalue Amax 
is not in the Krylov space (see Section E|- Furthermore, the accuracy of the esti- 
mation of the rightmost eigenvalues is only about 10"'^. If one wants to determine 
the bifurcation point more accurately one has to apply a more efficient algorithm. 
For instance, one may adapt the block Arnoldi method [78] to the Cayley-Arnoldi 
algorithm. In the following the developed algorithms are utilized to study important 
open questions related to the dewetting of a thin film (Section [U and the depinning 
of a pinned drop from a wettability defect (Section ^ . 

5. Short-time dynamics and coarsening for dewetting thin films. The 

dewetting process of a thin film can be initiated by two different mechanisms: Either 
via a surface instability (spinodal dewetting) or by heterogeneous nucleation at finite 
size defects [211 ESI [HI ED]. The only process that can initiate dewetting for metastable 
flat fllms is nucleation. Such a fllm is linearly stable and the flnite perturbation has 
to be larger than a critical threshold given by an unstable steady state solution. The 
situation is, however, different for linearly unstable fllms. In general, there exists a 
critical wavelength Ac = 27r/fcc. Any perturbation associated to a wavelength A > Ac 
grows exponentially in time with a growth rate (3 = —■m{ho)k'^{k'^ — k1). The re- 
sulting short-time dewetting structure consists of a regular drop or hole pattern. The 
characteristic distance between them equals normally the wave length Am = V^Ac cor- 
responding to the linear mode growing with the largest rate Pm- However, whether 
this instability mode represents the dominant mechanism depends on the character 
of the primary bifurcation: Deep inside the linearly unstable regime the bifurcation is 
supercritical and the fllm develops a surface instability with the properties described 
above. However, closer to the metastable region the primary bifurcation is subcriti- 
cal. The present threshold solutions offer as second pathway of evolution a nucleation 
process as in the metastable regime [1011 1103[ [89t H] . The latter process is normally 
dominant if the growth rate of the threshold solution is much larger than the linear 
growth rate Pm of the flat fllm. Under these circumstances, details of the dewetting 
process strongly depend on experimental conditions (amount of defects, roughness, 
noise). When nucleation dominates one expects larger drops or holes that are ran- 
domly distributed. We will use our algorithms to investigate (i) the dominance of 
either instability- or nucleation-triggered dewetting in the Hnear unstable thickness 
region in the 2d case; and (ii) quadratic and hexagonal steady state solutions in the 
2d case and the character of their primary bifurcations. Both questions were not yet 
studied in the literature. 

The short-time dewetting dynamics is followed by a very slow coarsening process. 
Indeed for a horizontal substrate, the system is variational and the only globally stable 
solution corresponds to the absolute minimum of the underlying Lyapunov functional. 
The corresponding solution is a single drop coexisting with the precursor fllm. The 
coarsening process corresponds to a cascade of two- (and three-) droplet mergers. This 
mass transfer can occur by two mechanisms based on a volume mode and a translation 
mode, respectively [IS]. In the volume mode all mass of one droplet flows through the 
precursor fllm into the other droplet. The centers of mass of the droplets do not move. 
In contrast, in the translation mode the entire droplets move, approach each other 
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and merge. Both modes are related to the Goldstone mode of translational invariance 
for a single liquid front on an homogeneous substrate 09]. Their stabilization by a 
substrate heterogeneity is discussed in Ref. [91]. For our case, it implies that during 
coarsening the main part of the dynamics occurs in the contact line region of the 
drops. Furthermore, because this motion is related to the translation invariance the 
corresponding eigenvalues are close to zero, explaining the slowness of the process. 
Note that, from the point of view of a dynamical system a coarsening step can be 
interpreted as a heteroclinic connection between steady states. 

5.1. The one-dimensional case. In order to compare our results for the one- 
dimensional case to the literature, simulations are performed at parameter values used 
in Ref. |103j employing disjoining pressure (I). The starting profile is a flat film with 
a small localized defect at the center. We utilize the Cayley transform with a regular 
mesh with = 1600 for a domain size I = 16Am (100 discretization points per Am)- 
The cases of dominating surface instability and dominating nucleation only differ by 
the initial film thickness {H = 2.4 and H ~ 3.2, respectively). 

Figures ??(a) and (b) give space-time contour plots of the long-time evolution 
in the two cases. In the case of a dominance of the surface instability (Fig. ??(a)), 
the initial growth phase {t < ISr™) clearly results in a regular droplet assembly 
corresponding to the fastest linear mode of wavelength Am- In particular, we find 16 
drops at equal distances. The profile and the evolution of the norm and the relative 
energy are in good agreement with the results in |103j (their figures 14(a) and 16(a)). 
For H = 3.2, the growth of a critical disturbance and the resulting hole is faster 
than the surface instability as expected. Further holes are subsequently nucleated 
by secondary nucleation events outside the rim around growing holes. The resulting 
short-time structure consists of larger drops than in the surface instability regime. In 
particular, before coarsening sets in we have only 9 drops. The duration of the initial 
'rupture' in ??(b) is in good agreement with Fig. 16(b) of Ref. |103j . 

With the here presented method one is able to study the long-time coarsening far 
beyond the results of, e.g., [1021 191]. Note the logarithmic time scale in Fig. ??. One 
can clearly see that for a drop pattern both above mentioned coarsening mechanisms 
play a role. These mechanisms are in agreement with the one discussed in Ref. |36) . 
The merging of drops does not occur continuously slow but an extremely slow start of 
the process is combined with a very fast culmination. In consequence, the evolution 
of the energy in time shown in Fig. ?? shows long plateaus connected by 'jumps' 
resulting from the fast part of the transition to the coarser state. The combination of 
slow and fast dynamics is well simulated using our adaptive time-step method. 

According to figure 11 the coarsening process is in the nucleation case slower than 
in the one dominated by the surface instability. Furthermore, it proceeds mainly via 
the mass transfer mode. Finally, two [five] drops remain in the surface instability 
[nucleation] case. In theory, coarsening should continue, but the evolution becomes 
so slow that we reach the limit of numerical accuracy, i.e., the eigenvalues related to 
the coarsening modes become smaller than the numerical accuracy. In particular, for 
leading eigenvalues smaller than 10"'', the numerical noise is not negligible and we 
are not able to observe the next coarsening step. 

In conclusion we can state that the developed algorithm well simulates the short- 
and long-time behavior in the Id case. In particular, the short-time evolution agree 
well results obtained using a semi- implicit scheme with a constant time-step |103) . 
In contrast, here the time-step varies by 6 orders of magnitude allowing us to study 
the long-time coarsening. We reach the numerical limit during late coarsening stages 
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because the leading eigenvalues are then too close to zero. 

5.2. The two-dimensional case. After utilizing the Id case to show the re- 
liability of our algorithm, we next employ the algorithms to study dewetting in the 
2d case. All linear stability results introduced above still apply. In particular, the 
parameter regions where films are linearly unstable remains the same as do the most 
dangerous wavelength A™ and corresponding growth rate However, in contrast 
to the Id case, 2d patterns involve two wavevectors ki and k2 that can lead to a 
number of different periodic steady state patterns [2^. Here, we track square and 
hexagonal patterns by imposing a periodic 30 x 30 square domain and a 30 x 30\/3 
rectangular domain, respectively. Choosing the mean film thickness H as control pa- 
rameter we obtain the bifurcation diagrams in Figs. ??(a) and ??(a) for the square 
and hexagonal patterns, respectively. Note that the finite domain size results in crit- 
ical film thicknesses that differ from the values for an infinite domain. In particular, 
one finds the upper limit at H2 = 2.95 instead of the value He = 3.63 expected for 
an infinite domain. Fig. ??(a) shows that for the square pattern at H2 a subcritical 
branch bifurcates from the trivial one. It becomes stable and turns back towards 
smaller thicknesses at a saddle-node bifurcation at Hsn — He- In consequence, for 
H2 < H < He the fiat film is metastable and can only dewet via a nucleation process 
that allows it to overcome the unstable threshold solution. All steady state profiles 
above H w 2.5 correspond to hole patterns (see Fig. ??(b)). Decreasing H further 
eventually the structure becomes a pattern of drops. The change corresponds to the 
rather steep change in the norm at H ~ 2.2. At the transition both, holes and drops, 
show a square shape and form a checkerboard pattern (rotated by 7r/4, see states 5*2 
in Fig. ??(b)). At smaller H one always finds drop patterns. The branch turns again 
and becomes unstable at another saddle-node bifurcation before it joins the trivial 
solution subcritically at H Hi. 

In contrast to the square pattern, the hexagonal pattern does not display such 
a transition between drop and hole structures along a single branch. Interestingly, 
one finds two distinct branches (Fig. ??(a)): One consists of hexagonal patterns of 
holes whereas the other one shows only hexagonal drop patterns (Fig. ??(b)). In 
consequence, at each of the two bifurcation points at Hi and H2 two branches bifurcate 
from the trivial fiat film solution. At the lower (upper) critical thickness the drop 
branch bifurcates subcritically (supercritically), whereas the hole branch bifurcates 
supercritically (subcritically). The respective subcritical branch gains stability at 
a respective saddle-node bifurcation. Remarkably, the supercritical branches do as 
well correspond to linearly unstable solutions. At some distance from the respective 
primary bifurcation point they gain stability through symmetry-breaking bifurcations. 

Although the bifurcation analysis for a small domain is very instructive it does 
not allow to predict in a direct way which dewetting mechanism dominates for larger 
system sizes. Next, we employ our time-stepping algorithm to study the dewetting of 
a film in a large square domain. The use of different initial perturbations will allow 
to discuss the prevalence of either dewetting by nucleation or by surface instability 
inside the linearly unstable thickness range. As initial condition we chose a defect 
that shall compete with the surface mode that emerges from the initial film rough- 
ness. The defect is radially symmetric with a profile h{r) that corresponds to the h{x) 
in the Id case. We consider film thicknesses H = 2.4 and H = 3.2 that lead in the 
Id case to surface-instability-dominated and nucleation-dominated processes, respec- 
tively. First, we do not include any initial roughness (Fig. ??) and find, in consequence, 
that for both thicknesses the short-time evolution conserves the radial symmetry. In 
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the instability dominated case (Fig. ??(a)) the distance between consecutive max- 
ima corresponds as expected to the fastest growing wavelength A™, whereas in the 
nucleation-dominated regime (Fig. ??(b)) the typical distance between rings is about 
twice Am. The radial symmetry is broken when the boundaries of the square box are 
felt. This influence kicks in earlier in the nucleation-dominated regime. Coarsening 
sets in after the short-time evolution. In the radially symmetric part it starts at the 
center and proceeds through a cascade of ring contractions. In the surface-instability- 
dominated case, driven by Laplace pressure the innermost ring contracts and forms 
a drop at the center (t=15.4 in Fig. ??(a)) that subsequently attracts the next rings 
(t=122, t= 259 and t=467 in Fig. ??(a)). In the outer non-radially-symmetric part 
coarsening also proceeds: in the surface instability [nucleation] regime it is dominated 
by droplet [hole] fusion via the translation mode. 

The good conservation of the initial radial symmetry shows that we have only a 
very small influence of (numerical) noise. This explains why normally in dewetting ex- 
periments with very thin fllms - that are much more prone to thermal noise and other 
tiny perturbations - no such regular structures are observed. However, recent related 
experiments with electrically destabilized thicker fllms show a radial ring structure 
when an inhomogeneous electrical fleld is applied in such a way that it plays the role 
of the flnite size central circular defect (see Figs. 2 and 4 of Ref. [19]). There noise is 
much less important as in our Fig. ??. The spirals observed in Ref. [79] might result 
from a similar process, however, with a different initial defect. 

We expect the relative importance of the processes driven by nucleation and 
surface instability to strongly depend on the noise in the system as a stronger noise will 
favor the surface instability. To determine its influence we perform several simulations 
using different surface roughnesses for the initial fllm. In particular, we add noise of 

0. 1% (Fig. ??) and 1% (Fig. ??) relative amplitude to the initial flat fllm with the 
centered defect. 

As expected, the radial structure starting from the central defect has less time to 
evolve the larger the initial noise because the noise 'accelerates' the isotropic spinodal 
process. For the small 0.1% noise (Fig. ??) the radially-symmetric structure is still 
appreciable quite some time even into the coarsening process. However, coarsening 
eventually 'washes out' any memory of the initial defects. In particular, for the 
surface-instability-dominated case the second of the evolving rings is still circular (at 
t — 4.92 in Fig. ??(a)), but already the depression outside that ring has lost its 
complete radial symmetry. Instead it is a ring-like assembly of holes (at t = 6.19). 
In the nucleation case, only the flrst ring is circular (at t = 4.29 in Fig. ??(b)). The 
depression beyond it, however, emerges already as a circular pattern of holes, i.e., as 
a typical secondary nucleation pattern (cf. [1]). In both cases the initial defect has no 
further influence on the structure. The remaining area is covered by typical spinodal 
structures (see, however, the discussion of Fig. ?? below). For the stronger initial 
noise of 1% (Fig. ??) the initial defect has less influence on the resulting structure, 

1. e., the central radial structure is dominant during a very short time only (till about 
t = 3 . . A) and is later homogenized through coarsening. 

For comparison, we flnally simulate the dewetting starting with a flat fllm without 
defect but with an initial noise of 1% amplitude. A qualitative difference between 
H = 2.4 (Fig. ??(a)) and H ~ 3.2 (Fig. ??(b)) clearly shows already in the short-time 
behavior. At t/to = 3.8 some holes are already formed for the thicker fllm while for 
the thinner one only slight surface modulations appear. At larger times labyrinthine 
and hole structures evolve and coarse for H = 2.4 and H = 3.2, respectively. Late 
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in the process, the labyrinthine pattern transforms into a drop pattern while for 
H = 3.2 the holes remain. However, also in the latter case in a larger domain and 
at later times drops will prevail. Note finally, that the location of the larger holes 
for the respective second snapshot for the two time evolutions in figure ?? are almost 
identical. This proves that the numerical noise is negligible as we used identical 
initial noise in the two simulations. The interpretation of the described behavior is 
subtle: The different morphologies observed for H = 2.4 and H = 3.2 is not related 
to being in the surface-instability- or nucleation-dominated regime as there exists no 
initial finite disturbances. These are normal non- linear structures evolving through 
a surface instability, i.e., changing from drops to labyrinths and finally to holes as 
one increases the film thickness (cf. Refs. [II1II2])- However, the higher speed of the 
non-linear process at if = 3.2 might well be related to the fact that we are in the 
nucleation-dominated regime, i.e., the surface instability allows modulations to grow 
that then act as nucleation sites as the faster nucleation modes kick in. 

6. Depinning of a drop on a heterogeneous substrate . The second prob- 
lem we focus on is the depinning of a pinned droplet. If one applies a lateral force on 
film/drop on a homogeneous substrate [P 7^ in Eq. I|2.ip ]. one only finds traveling 
surface waves or sliding drops [1021 [9i] . No steady-state solutions exist, except the 
trivial fiat film. On a heterogeneous substrate the situation is different because the 
heterogeneity (e.g., chemical or topographical defect) can pin a drop or a hole. Wet- 
tability defects are modeled by varying the disjoining pressure along the substrate. 
Here, we consider heterogeneities in the form of stripes, i.e., we modulate the disjoin- 
ing pressure only in one spatial direction (here the x-direction) . We choose the defect 
profile i{x) based on Jacobi elliptic functions as described in [971 [96]. The specific 
used wettability profile is represented in the lower part of the panels of Fig. ??. The 
used disjoining pressure (HI) is given by Eq. (|2.4p with b = 0.1. The parameter e 
represents the amplitude of the heterogeneity, i.e., the wettability contrast. If e > 
[e < 0] the defect is less [more] wettable than the surrounding substrate, i.e., the 
defect is hydrophobic [hydrophilic] . 

When the lateral driving, e.g., the inclination of the substrate, is increased the 
pinned droplets are deformed and their center of mass slightly shifts until at a critical 
driving Pc the droplet depins. The analysis of the steady states and the depinning 
bifurcation is tackled using the continuation approach developed above in Section IH 
In the Id case, results are already available and can be found in Ref. [971 ES] were 
the continuation package AUTO [30] and an explicit time-stepping algorithm with 
adaptive time-step were used. We employ this case in Section 16.11 to validate our 
continuation code. As the results for identical parameters are entirely identical to 
the ones of Ref. |97[ 196). we here present results for larger drops. In the subsequent 
Section 16. 2p we explore the physically more realistic two-dimensional case which can 
not be treated using AUTO because the governing equations are not equivalent to an 
ODE system. 

6.1. The one-dimensional case. First we consider shortly the Id case. The 
domain size and the volume of liquid volume are fixed at L = 50 and V = To, respec- 
tively (i.e., mean film thickness is if = 1.5). For a homogeneous substrate without 
lateral driving (e = and P = 0), the critical wavelength for spinodal dewetting is 
Ac — 15, i.e., for the chosen domain size there exist at least steady states containing 
one, two or three drops. They bifurcate from the trivial fiat film solution at nAc with 
n = 1, 2, 3, respectively. If the primary bifurcation is subcritical there might be more 
solutions (cf. Ref. [95j). 
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To determine the various steady-state solutions for a heterogeneous substrate 
without lateral driving we use continuation when varying the heterogeneity strength 
e for a heterogeneity period equal to the domain size. The initial solutions are readily 
found starting from a flat fllm solution without heterogeneity (Fig. ??, cf. as well 
Ref. [HI])- Branches of steady-state solutions cross the axis e = three times, corre- 
sponding to the one-, two- and three-drop solution for a homogeneous case that we 
introduced above. For a strong heterogeneity (|e| sufficiently large) only the steady 
state remains that corresponds to a single drop. This is the most interesting solution 
for a study of depinning. 

However, for smaller |e| up to seven steady states can exist corresponding to var- 
ious stable and unstable one-, two, and three-drop equilibria. As Ref. [97] studied a 
smaller domain (L = 25) their Fig. 6 shows only one crossing of the axis e = 0. Study- 
ing the stability of the equilibria one flnds that only the uppermost branch corresponds 
to stable solutions. All other branches terminate in saddle-node bifurcations. 

Next, all stable and unstable steady-state branches are tracked when increasing 
the lateral driving force P. Bifurcation diagrams and selected proflles of pinned 
droplets are given in Figs. ?? and ??, respectively. The various branches found for 
P = continue to exist for small driving as 'pinned solutions'. However, at critical 
values of the driving most cease to exist as they annihilate each other. Physically 
speaking, the heterogeneity can not retain the drops any more and they start to sHde, 
i.e., they depin. Beyond a critical value Pc no steady drop does exist as even the 
stable single drop depins (Fig. ??). The bifurcations where solutions annihilate are 
for the shown cases SNIPER (Saddle Node Inflnite PERiod) bifurcations. At Pc a 
branch of space- and time-periodic solutions emerges (not shown) that corresponds to 
sHding drops that are modulated by the heterogeneity profile. The temporal period 
diverges as one approaches Pc from above and the drop motion resembles stick-slip 
behavior (cf. Refs. [971196]. 

Note, however, that selected branches exist for all P. We can distinguish two 
cases: (i) small amplitude branches showing no bifurcation start at P = as the 
unstable three-drop state, e.g. the lowest curve for e = ±0.3 in Fig. ??; and (ii) the 
branch of the initially (at P = 0) stable one-drop solution which undergoes two or 
more saddle-nodes, e.g., the curve for e = ±0.7 in Fig. ??. In both cases for P > Pc 
the branch consists of small amplitude solutions, i.e., it corresponds to film fiow where 
the steady film surface is modulated by the heterogeneity. We find all such solutions 
to be unstable. 

The results obtained up to here indicate that our continuation algorithm is well 
capable to follow stable and unstable steady states in the Id case. The saddle-node 
bifurcation has been well detected and as expected no numerical stability problems 
have been encountered near the bifurcation. Next, the continuation algorithm has to 
prove its capabilities in the 2d case. 

6.2. The two-dimensional case. In the 2d case we consider the depinning of 
drops from stripe-like wettability defects. In particular, we employ an x-dependent 
heterogeneity profile and choose the lateral forcing to be as well in x-direction. In 
this way a hydrophobic stripe blocks a drop at its front end whereas a hydrophilic 
one holds it at the back end. The steady-state droplets are continued for increasing 
lateral driving force P. Fig. ?? presents the bifurcation diagram for a single stable 
drop blocked by a hydrophobic stripe with e = 0.3. Stable (unstable) steady states are 
given as solid (dashed) lines: The stable drop blocked by the defect increases its norm 
with increasing driving as it is 'pressed' against the defect. The drop finally depins 
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at a saddle-node bifurcation {Pc ~ 0.003) where the stable branch turns back and 
loses stability. Time simulations beyond depinning show a typical stick-slip motion 
of drops with a period that diverges when the bifurcation is approached from above. 
The time averaged norm for selected values of P is indicted by triangles in Fig. ??. 
The results indicate that depinning occurs via a SNIPER bifurcation. 

Examples of steady stable and unstable drop profiles are given in Figs. ??(a) and 
(b), respectively. The stable drop sits behind the Hne of minimal wettability of the 
defect. As it is squeezed against the defect by the driving force it has an oval shape. In 
contrast, the drop that represents an unstable equilibrium, has a 'forward protrusion' 
that crosses over the minimum of wettability. This solution represents a threshold 
for depinning for P < Pc, i.e., a small perturbation will either let the drop retract its 
advancing protrusion to again settle behind the defect or trigger a depinning event 
that sends the drop sliding down to the next defect where it is pinned again. 

Note, that at large values of P, the norm of the remaining solutions on the lower 
branch does not tend to zero as in the Id case where one finds an (unstable) modulated, 
nearly fiat film. Indeed, the solution becomes a spatially modulated steady rivulet: a 
solution with small variations in the x-direction, with profiles in the y-direction that 
are similar to Id drop solutions. 

Let us, finally, come back to the computational problem. Apart from the increase 
of the system size N, in the 2d case one encounters a new difficulty related to the 
translation symmetry in y-direction. It implies that for each solution there exist 
a periodic orbit of solutions obtained by translation in y-direction that has to be 
avoided by the continuation algorithm. Neglecting numerical noise, the solutions 
possess a left-right refiection symmetry y —y (Fig. ??). Therefore, the Jacobian is 
an equivariant operator for this left-right refiection. Thus, the action of the Jacobian 
on left-right symmetric vectors results in vectors with the same symmetry. This 
eflFectively excludes any translation in the y-direction. However, even small numerical 
noise becomes relevant when the solution is close to the trivial fiat film state and the 
leading eigenvalues are very small as well. Indeed, that is exactly what we observe 
numerically in the above example. For P > 0.1 the solutions correspond to almost 
fiat rivulet-type solutions and the continuation algorithm might stay on the orbit of 
solutions related by translation. This results in a mere translation of the rivulet in 
the y direction. The problem can be easily overcome by fixing the maximum at a 
particular point. However, this might cause problems in more complex cases when 
various maxima may coexist. To avoid any ambiguity, we prefer to use the same 
technique as for the problem of mass conservation. In particular, when applying 
the Krylov procedure, we project the basis vectors orthogonally to the translation 
mode dyho. This does not change the structure of the algorithm and is performed at 
negligible cost (see Section [A|. 

7. Conclusion. We have developed (i) a time integration scheme based on expo- 
nential propagation and (ii) a continuation algorithm employing the Cayley transform 
for the non-Hnear lubrication equations. These equations contain difi^erential opera- 
tors till fourth order. To avoid severe stability restrictions on the time-step, the linear 
term may be treated implicitly. However, for the lubrication equations the difiiculty 
is to find a relevant linear operator. In consequence we use the Jacobian matrix at 
each time-step. In this framework, an exponential propagation scheme is more effi- 
cient than a semi-impHcit scheme [inil04| . The method is based on an exact solution 
of the linear problem for each time-step. This involves the determination of the ex- 
ponential of a matrix. To do this in an economic way, the linear operators may be 
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reduced to small Krylov subspaces. However, we find that for the lubrication equation 
(that contains the fourth order Bilaplacian operator) the necessary dimension of the 
Krylov subspaces is typically about 100 - much larger than used in Ref. [34] for a 
problem involving only a second order operator. 

To accelerate the convergence, we have coupled the Arnoldi algorithm with the 
Cayley transform that is performed using an ILU-factorization. This Cayley-Arnoldi 
method accelerates the convergence of the Krylov reduction and, in consequence, 
allows to take larger time-steps. Furthermore, the leading eigenvalue is well estimated 
and can be used to effectively adapt the time-step to the changing characteristic time 
scale of the dynamics. This allows for a major improvement in simulations of the very 
slow coarsening dynamics of one- or two-dimensional systems [1031 [T2l [9T| [84] . 

We have also developed an algorithm for the continuation (or path-following) of 
steady states that is based on the tangent predictor - secant corrector scheme. Both 
tasks - time stepping and continuation - can be performed using the same Cayley- 
Arnoldi algorithm. The advantage of this approach is the possibility to perform all 
the bifurcation analysis tasks at the same time. This includes the computation of the 
kernel of the .Jacobian for the bifurcation detection, and stability analysis (evaluation 
of leading eigenvalues) . 

The developed algorithms have been used to study the bifurcational structure 
and time evolution of (i) dewetting thin films and (ii) droplets depinning from a 
heterogeneity for physically two- and three-dimensional settings, that correspond to 
Id and 2d equations, respectively. 

For the dewetting film we have investigated different pathways for the initial 'rup- 
ture', i.e., the short-time evolution. The focus has been on the competition of the 
surface instability and the nucleation at defects. The long-time coarsening dynamics 
has as well been studied. For the Id case, we could follow the coarsening process 
till lO^Tm where Tm is the characteristic time of the surface instability. In the two- 
dimensional case, we have found that the short-time dewetting process induced by 
a circular defect conserves the radial symmetry until disturbed by the boundaries 
of the square domain. This indicates the very small influence of round-off errors in 
our algorithm. In this case the long-time coarsening proceeds via a cascade of ring 
contractions. However due to thermal fluctuations such regular structures are not ob- 
served in dewetting experiments with very thin fllms. However, they are observed for 
dielectric fllms in a capacitor when an inhomogeneous electrical fleld is applied in such 
a way that it plays the role of the flnite size central circular defect [1^ . Adding noise 
to the initial conditions we recover the 'classical' labyrinthine dewetting structures 
deep inside the linearly unstable region or hole patterns for thicker fllms. Using the 
Cayley-Arnoldi method for a system size of about 10 Am x lOAm, we have been able to 
simulate the coarsening dynamics till reaching a single drop, i.e., the globally stable 
solution (Fig. ??). Beside the time evolution we have employed continuation to study 
for the flrst time 2d steady state solutions corresponding to square and hexagonal 
arrangements of drops or holes. 

Furthermore, we have applied our framework to study the depinning of ridges 
(Id case) and drops (2d case) that are pinned at substrate heterogeneities. The 
pinned drop solution has been followed in various parameters using the developed 
path-following algorithm. In particular, we have used the wettability contrast related 
to the amplitude of the proflle of the heterogeneity e and the lateral driving force P. 
The one-dimensional case we have used as a benchmark for our algorithm and have 
successfully reproduced results obtained in Refs. |97[[96] using the package AUTO |30j . 
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We have as well been able to track two-dimensional stable and unstable steady drop 
and rivulet solutions. This has been supplemented by a study of the time-dependent 
solutions beyond depinning. This has led us to the conclusion that depinning occurs 
via a SNIPER bifurcation and that beyond (but close to) the bifurcation the droplets 
move showing a stick-slip behavior. 

Although our approach improves the time integration and path-following for lu- 
brication equations, the constant equidistant finite difference spatial discretization re- 
mains very basic. The weakness of such a regular discretization appears, for instance, 
when tracking large drops pinned at hydrophobic defects. For increasing driving force 
the stable drop becomes strongly localized at the defect. However, the unstable solu- 
tions for the same driving force is very close to the stable one (as already observed in 
the 2D case for very large drops - see Fig. 29 of Ref. [97]). To clearly distinguish the 
two solutions numerically a higher accuracy is required. This can be achieved using 
an adaptive mesh. In our particular system, for instance, the continuation procedure 
applied to a drop on a 50 x 50 square domain breaks down near the depinning bifur- 
cation. This accuracy problem is similar to the one encountered in the last stages of 
coarsening in Section ISTTl where a mesh refinement is required near the droplet edges. 
An adaptive mesh would not only decrease the size of the discretization space but we 
expect that it would also improves the conditioning of the Jacobian matrix. 

The physical situations we have presented here have been restricted to the specific 
lubrication equation for films and drops on partially wettable homogeneous and het- 
erogeneous substrates. Beside the curvature pressure the disjoining pressure has been 
the only term resulting from the underlying free energy. However, the method can be 
applied to various thin film systems involving other contributions to the free energy. 
This includes, for instance, thin films of dielectric liquid in a capacitor [56 ^ \6U \ flOS) . 
heated thin films unstable due to a long- wave Marangoni instability |67[ 1106^ \T2\ [95] , 
films with an effective thickness dependent surface tension caused by high-frequency 
vibration [53l 1100) and film flows in a time-dependent ratchet potential [171 HB] . 

The algorithms can also be applied to multilayers problems in the long-wave 
framework like, for instance, in bilayer dewetting [TTJ [33l [21 [72] or to systems involving 
a single thin fllm equation coupled to a (reaction-) diffusion equation for a reactive 
[23l [68] or non-reactive |44[ [58] surfactant or adsorbate at the substrate |93[ [45] . For 
those problems the overall structure of the equations does not change. Only the system 
size N is multiplied by the number of equations. The time-stepping and continuation 
code can be applied without major change. 

In general, we expect the method to be relevant as well for closely related evolu- 
tion equations like, for instance, the (driven) Cahn-Hilliard equation |17[ [37] or the 
Kuramoto-Sivashinsky [5H[T5] equation which both contain the Bilaplacian operator. 
However, in contrast to the lubrication equation for those equations the mobility func- 
tion m{h) multiplying the Bilaplacian is a constant (although this is an approximation 
in the case of the Cahn Hilliard equation). As mentioned at the begin of Section [3l 
the presence of the non-constant mobility function has been one reason for our choice 
of the Jacobian matrix as Hnear operator at each time-step. This implies that for the 
constant-mobility Cahn-Hilliard and the Kuramoto-Sivashinsky equation, this choice 
is not crucial and a classical semi-implicit scheme with L = as Hnear operator 
might result in a viable scheme. The efficiency of both time integration schemes has 
to be compared to decide which is the most powerful method. 

The presented continuation algorithm for steady-state solutions, can be improved 
in a straight forward manner by adapting it for stationary states, i.e., for travelling 
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waves or sliding drops. These can be seen as steady-state solution in a co- moving 
frame. For the case of a driving force in x-direction, solutions are also invariant 
w.r.t. translation in a;-direction. The resulting problems can be overcome using the 
same technique as above for the translational invariance in the y-direction. Thus, 
the Cayley-Arnoldi method can still be applied without major change. The extended 
algorithm would be applicable to the study of sliding drops on inclined homogeneous 
substrates. In particular, the onset of the morphological transition between single sta- 
ble sliding drops to drops that emit small droplets at the back would be a worthwhile 
subject for further study [69l [86] 

Appendix A. Krylov reductions. 

A.l. Arnoldi-Krylov. The approximation of Vg = G(Jt)6 and Ve = exp(jT)c 
constitutes a crucial step of our time integration algorithm (end of Section [3]). Here, 
we propose to use a Krylov reduction as it is usually employed for sparse operators. 
One aims at obtaining an accurate approximation so that the time-step is only limited 
by the order of the scheme. Because the approximation of both vectors-Wg and ■^e^can 
be performed using the same technique, we only focus on Vg (corresponding to the 
second order linear scheme Eq. p.l2p ). The main idea of the Krylov reduction is that 
the series of subspaces 

(A.l) = span{6,J6,j26,..., J™-ifo} 

converge to a finite dimension Krylov subspace Km which contains Vg. Obviously, 
it is only an efficient method if the subspace Km is a good approximation of Vg for 
m <^ N . The Arnoldi method is a standard method to construct an orthonormal 
basis Vm of the subspace K^.- The approximated Jacobian matrix, denoted J™, is 
then a TO X m upper Hessenberg matrix: 

(A.2) J™ - V^3Vm. 

Relation l|A.2p can be used to approximate 

(A.3) Vg = G(Jt)6 ~ y„,G(J,„r)y„\6. 

Since the dimension m is small, the classical QR algorithm is a reliable and efficient 
method to diagonalize J„i to obtain the matrix Dm and then to compute G(3mT). 
The sought for approximation of the vector Vg is given by 

(A.4) Vg ~ Vm'P^G{DmT)V^V^b 

where Pm is the matrix of the eigenvectors of Jm • The columns of the rectangular 
matrix (KnPm) are the Ritz vectors, i.e., the approximated eigenvectors of J. Since 
the dynamics is dominated by the rightmost eigenvalues of J a good Krylov approxi- 
mation impHes that the Ritz vectors are close to the rightmost eigen-directions. We 
call them the 'wanted' eigendirections. However, it is known that the Krylov- Arnoldi 
algorithm converges first to the eigenvalues of largest modulus. The latter are situated 
in the leftmost spectrum and are the main reason for stability problems. We call them 
'unwanted'. Our aim is to improve the Krylov- Arnoldi method by adapting algorithms 
used to estimate the rightmost spectrum. To do so we transform the spectrum of J 
in such a way that the wanted eigendirections are associated to the eigenvalues of 
largest modulus. Thus, by applying the classical Krylov- Arnoldi to this transformed 
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operator the wanted eigendirections are selected after a few steps. There exist two 
main transformations for asymmetric sparse systems — Chebyschev acceleration and 
shift invert Cay ley transform — whose efHciency will be discussed next. 

Note, that the Krylov reduction allows to introduce supplementary requirements 
on the discretized space in a simple way. For instance, in order to preserve the 
volume during appHcation of the Arnoldi algorithm the direction corresponding to a 
variation of volume is suppressed. Therefore, one ensures that the height h remains 
in the Euclidian space {H + Eq} (Section 12. 2p during the time-stepping. In a similar 
manner, the directions related to translation invariance may be suppressed during the 
continuation algorithm to avoid problem as mentioned at the end of Section 16.21 

A. 2. Chebyschev acceleration. The approximation (Eq. I|A.3P ) can be inter- 
preted as a polynomial approximation since the basis Vm is a sum of elements of Km 
(Eq. {AH)) [77]: 

(A. 5) VG^Pm{JT)b 

where p„i is a polynomial of degree to — 1. The principle of the Chebyshev acceleration 
is to seek an optimal polynomial to accelerate the convergence. In [ST] it is shown that 
scaled and translated Chebyschev polynomials have certain optimal convergence prop- 
erties. The main reason of the success of the Chebyschev polynomial is the possibility 
to decrease some 'unwanted' eigendirections contained in an ellipse. This property is 
employed to compute the rightmost eigenvalues of large sparse non-symmetric matri- 
ces [39t 178). However, in our case the application of the algorithms presented in [59] 
does not significantly improve the convergence of the Krylov approximation. Indeed, 
according to [78] the rapidity of the convergence is directly affected by the accumula- 
tion of the rightmost eigenvalues. In consequence, a polynomial approximation does 
not result in the improvement of the Krylov reduction. 

A. 3. Cayley transform. Unlike the Chebyschev method, the Cayley transform 
is not polynomial but rational. Let us introduce the transformed operator C: 

(A.6) C = J-i = (J-cI)-i 

where c is a real constant that has to be chosen. The matrix C contains the eigen- 
vectors of J but has a different spectrum. If c is chosen to be larger than the leading 
eigenvalue Amax of J, the spectrum of C falls into the band [(Re(Ainax) — c)~^;0[. 
In consequence, the wanted eigen-directions of J correspond to the eigenvalues of 
C with the largest modulus. Therefore, the orthonormal basis Vm constructed us- 
ing the Arnoldi procedure within C should converge after a few steps to the wanted 
eigen-directions. One introduces the approximate operator Cm 

(A.7) Cm = V,iCVm 

and diagonalizes Cm using the QR method (similar to the Arnoldi-Krylov reduction): 

(A.8) Cm Pm-D cayley -Pm * 

Finally according to the definition IjA.GP and the Krylov reduction (|A.7|) we obtain 
the approximation of Vg 

(A.9) Vg = G(Jr)& c K.PmG(D;^^,^y + cI)P^^V:;,b. 
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The efficiency of the algorithm depends strongly on the choice of c. Indeed if c » 
Re(Aniax), then the Cayley transform is not relevant any more since C ~ — c^^I. In 
addition, if c is close to an eigenvalue of J, the operator J— cl becomes singular and the 
method diverges. Numerical calculations show on the one hand, that for moderately 
large values of c (as compared to Re(Aniax)) the accuracy decreases notably. On the 
other hand, a choice of c close to Re(Ainax) leads to very accurate results even if 
c < Rc(Aniax)- In consequence, a good choice is to take a constant c'^ at time-step 
k that is slightly larger than Re(A^~;^) estimated in the previous time-step fc — 1. 
The constant c'^ could then be smaller or larger than Re(A^jjx). In rare cases it might 
be very close to an eigenvalue of the Jacobian matrix J. However, such a degeneracy 
would be detected automatically by the presence of huge rightmost eigenvalues of the 
matrix C. So, in that singular case the time-step is performed again with a larger 
value of c. 

The difficulty of the method remains the evaluation of the action of C, de- 
fined as the inverse of Jc, on the vector b. It is obtained utilizing an Incomplete 
LU-factorization (denoted ILU) on Jc which is a powerful method for sparse band- 
matrices. The cost of this factorization is about 0{N^/'^). All other operations are 
0{N) we expect that the ILU slows down the algorithm for large systems. 

A. 4. Comparison. We next compare the Krylov-Arnoldi and Cayley- Arnoldi 
methods for the estimation the vector Vg. Even though the Krylov reduction does 
not depend on the chosen time-step r, we expect that the approximation of Vg does. 
We aim at a Krylov reduction that is accurate enough to not add a restriction on the 
time-step t. For a second order linear scheme Eq. I|3.12p the relevant time scale is 
about the inverse of the leading eigenvalue, denoted t\. This value is employed in the 
different numerical convergence tests. 

In the one-dimensional case, Vg may be computed by a direct method as, for 
instance, a QR-diagonalization. The latter is taken as reference value Wrof for the 
relative error. Unfortunately, for the two-dimensional case, the system may be large, 
i.e. N = O(IO^), and a QR-diagonalization dramatically increases the CPU cost and 
memory requirements (being proportional to A^'^). Thus, in this case, the reference 
solution Wrcf is computed using the Cayley- Arnoldi method for a Krylov subspace 
dimension m large enough to obtain convergence. 

Figs. ?? and ?? present selected results for the convergence of the time-stepping 
scheme for the dewetting problem in the Id and 2d case, respectively. Both panels 
(a) shown an extremely slow convergence of the classical Krylov method. Krylov 
subspaces with dimensions m w 100 have still a relative error of 10~^. An accuracy 
of about 10~^ is obtained by a 800-dimensional Krylov subspace (Fig. (??)(a)). In 
contrast, with the Cayley-method, the same accuracy is obtained with Krylov sub- 
spaces of much small dimensions m ^ 10 . . . 30. For a time-step r smaller than the 
characteristic time r^, the convergence is notably improved for the classical Krylov 
algorithm only. This indicates that for this Krylov reduction the eigenvalues are badly 
approximated contrary to a Krylov algorithm using the Cayley transform. We find 
that the efficiency of the two methods is equivalent for N ~ 10^. By 'efficiency' we 
mean the ratio between the time-step r and the CPU time cost for a relative error of 
10-6. 

In conclusion, the Cayley transform emerges as a powerful method to perform an 
accurate Krylov reduction. However, its efficiency is limited by the ILU-factorization 
which considerably slows down the time-step for large system N > 10^. In this case the 
simple Arnoldi procedure is preferable and the accuracy of the time integration scheme 
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is determined by the Krylov approximation. Finally, the Cayley-Arnoldi reduction 
may be applied to higher order schemes such as those presented in [HI 1104) . Since it 
is not required to use another ILU-factorization this can be done at negligible cost. 

Appendix B. Discussion of algorithms. 

B.l. Algorithms. As the time-stepping algorithm based on the simple Arnoldi 
procedure is essentially as in [31], we only present the one that employs the Cayley 
transform. Since, the Cayley-Arnoldi method allows to estimate the leading eigenvalue 
Am which can be used to determine a time-step Tk independently of the previous time- 
step Tfc_i. The relative profile variation defined by 



can be approximated using Eq. I|3.12p with the Jacobian matrix J replaced by the 
scalar Am- Then for a given variation Cvar we obtain the time-step 



In the present work we fix the typical value of evar at about a few percent. The relative 
error er resulting from the estimate (|B.2p is rarely larger than a 10~^. Incorporating 
such a procedure to adapt the time step we obtain the computational schemes: 

Time-stepping with Cayley-Arnoldi method 
. We impose 

• the computational grid {N, 6x) 

• and the tolerance of the relative error etoi ', 
Start with ho,tQ] and 

1. Determine b — F{ho, x),J~ DhF{hQ)\ 

2. Compute the matrices L, U of the ILU-factorization; 

3. Construct Krylov subspace associated with h and J, determine the matrices 
K»i,Pm,Dm using the Cayley-Arnoldi method; 

4. Determine the leading eigenvalue A™ as good estimate, then deduce the time- 
step using Eq. (|B.2P : 

5. Compute u^^^ given by Eq. I|3.12p : 

6. Non-linear correction: determine the matrices Vm-, Pm, Dm using the Cayley- 
Arnoldi method; 

7. Evaluate the relative error e using Eq. p.lOp . distinguish two cases: 

(a) If < etoi the time-step is performed, 

(b) Otherwise decrease the time-step r t 12 and restart at step (5). 
Path- following with Cayley-Arnoldi method 

. We impose 

• the computational grid (iV, 5x) 

• and the tolerance of the relative error etoi\ 
Start with ho,po; and 

1. Determine b = -DpF{ho,po), 3 = DhF{ho,po); 

2. Compute the matrices L,U oi the ILU-factorization of J; 

3. Construct Krylov subspace associated with b and J, determine the matrices 
Kn,Pm,Di„ using the Cayley-Arnoldi method; 

4. Determine the tangent direction {ut, dpt), let hi — Hq+Ui and pi — po +dpt; 

5. While \\F{hk,Pk,x)\ \ > etoi we 



(B.l) 




\M 



(B.2) 
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(a) Determine c = ~F{hk,Pk,x), d = ~DpF{hk,Pk, x) and J = DhF{hk,Pk, x); 

(b) Compute V,n, Pm, using the Cayley-Arnoldi method (with L, U ma- 
trices of step (2)); 

(c) Inverse the equations l|4.15p and (|4.16p with condition (|4.18p . 

B.2. Comparison to 'classical' algorithms. The emphasis of both algorithms 
Ues on the Cayley-transformation of J using the ILU factorization. As the complexity 
of this task is 0(A^3/2-) Qp-^ ^^^^ 

is a more important issue as for other schemes in 
the literature, in particular, the time-stepping schemes. We may ask us how relevant 
is our choice? 

If one wanted to implement a classical implicit time integrator one would need to 
solve a linear system involving the matrix J. For instance, a backward Euler scheme 
leads to the system {I + 3t)u = br, where 6 is a known vector. A typical choice for the 
latter is the use of an iterative Krylov method as, e.g., the GEMRES method which 
is a good candidate for unsymmetric matrices. The efficiency of this method does 
strongly depend on the knowledge of an effective preconditioner. Without precon- 
ditioner, the inversion for a simple semi-implicit scheme (backward Euler) converges 
very slowly and almost fails to obtain the wanted tolerance. As noted at the end 
of Section [2] no general preconditioner exists for J. Therefore an ILU factorization 
of the matrix appears to be the only systematic way to construct an effective pre- 
conditioner. Therefore, a clear analogy exists between the exponential scheme and a 
semi-implicit scheme. On the one hand the Arnoldi algorithm to compute G{Jt) is 
equivalent to an iterative method without preconditioner to solve the linear system 
in an implicit scheme. On the other hand, the Cayley-Arnoldi algorithm to compute 
G(Jt) is equivalent to an iterative method with LU preconditioner to solve the linear 
system in an implicit scheme. Judging the complexity of the algorithm, a classical 
implicit and an exponential propagations algorithm are roughly equivalent schemes 
at the same order. However, exponential schemes seem to have two advantages: 

• The G and exp functions have a better leftmost spectrum filtering property 
than rational functions. 

• Literature results indicate that the Krylov techniques converge faster when 
employed for the evaluation of G(Jt)6 and exp(Jr) than when employed for 
the solution of a linear system (41] . 

Furthermore, does the particular Cayley-Arnoldi exponential propagation method 
used have the ability to give an estimation of the leading eigenvalues what the implicit 
method is not able to do. This information allows for an adaptive time-step that 
constitutes a major advantage when studying dynamics characterized by different 
time scales. In addition, it is an important information that allows to perform non- 
standard stability analysis as, e.g., employed by Miinch [63] to detect the formation 
of fingers during the dynamic opening of a single hole in dewetting. 

Although the shift Cayley-transform is a standard method to find the rightmost 
eigenvalues of a operator (see, e.g., Ref. [llOj ) its application in a continuation algo- 
rithm is less common. However, when solving a linear system within the Newton algo- 
rithm one faces the same problems as discussed above in the context of time-stepping. 
Thus, an effective solution method requires the use of the LU preconditioner. How- 
ever, our algorithm does not use the ILU factorization for a direct inversion, but rather 
to find the rightmost eigenvalues. The advantage of this approach is the ability to 
determine the tangent direction even close to a saddle-node bifurcation. Furthermore, 
at a bifurcation point, the directions of (different) bifurcating branches can be found 
during step 4 of the algorithm. 
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We emphasize that the numerical difficulties encountered in the time-stepping and 
the continuation task are both overcome using the Cayley transform. This is reminis- 
cent to the analysis of Tuckerman and Barkley |105j who point out that the classical 
implicit time-stepping scheme can be adapted to perform the bifurcation tasks. The 
developed algorithms overcome two main problems: (i) the operator has no 'simple' 
relevant linear part (different spatial scaling), and (ii) the "very bad" conditioning of 
the bilaplacian. For the current computers, the Cayley- Arnoldi method is the most 
efficient method for a moderate system size of = O(IO^). In consequence, the 
scheme is difficult to adapt it for three-dimensional PDE's. 
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